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Dissipative particle dynamics (DPD) does not conserve energy and this precludes its use in 
the study of thermal processes in complex fluids. We present here a generalization of DPD that 
incorporates an internal energy and a temperature variable for each particle. The dissipation induced 
by the dissipative forces between particles is invested in raising the internal energy of the particles. 
Thermal conduction occurs by means of (inverse) temperature differences. The model can be viewed 
as a simplified solver of the fluctuating hydrodynamic equations and opens up the possibility of 
studying thermal processes in complex fluids with a mesoscopic simulation technique. 



Dissipative particle dynamics constitutes a valuable tool for mesoscopic simulations of complex fluids. It was 
introduced by Hoogerbrugge and Koelman ffl , [Q and has received growing interest in view of its potentiality in the 
study of complex flow problems as those arising in porous flow 0, colloidal suspensions Q), polymer suspensions 
||, or multicomponent flows ||. The technique has received a considerable theoretical backup |t]]- (ll]] that provides 
solid ground for its use. 

DPD faces, however, a conceptual problem in that energy is not conserved: the dissipative particles interact with 
dissipative forces that depend on the relative velocities between particles. For this reason, a DPD system cannot 
sustain a temperature gradient as it has been pointed out by Marsh et al. JTo{ ] . Nevertheless, in the picture where the 
dissipative particles are understood as droplets or mesoscopic clusters of atoms fl|jL2|| it is apparent that the dissipated 
energy due to friction must be invested into increasing the internal energy of the cluster. Therefore, we propose in 
this letter to introduce an additional variable ti which is interpreted as the internal energy of each particle. Along 
this new variable we introduce an entropy variable Si — s(ei) and a temperature Ti = [dsi/dei]^ 1 . 

The variation of the internal energy must involve two different processes. On one hand, temperature differences 
between particles produce variations in the internal energies through "heat conduction" . On the other, the friction 
forces dissipate energy which is transformed into internal energy through "viscous heating" . Let us analyze each 
process separately. 



Let us assume for a while that the N dissipative particles of the system are at rest at positions r,. We formulate 
the following equation of motion for the internal energy of each particle 



This is a discrete fluctuating Fourier equation of heat conduction with k playing the role of a thermal conductivity. We 
note that the thermodynamic quantity conjugated to the internal energy is the inverse of the temperature rather than 
the temperature itself, = | — r j- 1 is the distance between particles i,j and the weight function to(r) determines the 
range of influence between particles. The random process cjij mimics the random heat flux that occurs spontaneously 
due to thermal fluctuations [juS| . This random "heat flux" is assumed to be antisymmetric under particle interchange, 
this is, — —qji in such a way that the total internal energy of the system is conserved 4r ti = if only a heat 
conduction process occurs. 

The Langevin equations (|l|) have a mathematical well-defined correspondence in the form of a stochastic differential 
equation (SDE), 
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I. INTRODUCTION 



II. CONDUCTION 
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where we have expressed the random "heat flux" q"ij from particle i to j in terms of the increments of the Wienner 
process dW^ . The function Co(rij ) provides the range of interaction and a is an overall noise amplitude. The increments 
of the Wienner process are antisymmetric under particle interchange dWfj = —dWL and satisfy the mnemotechnical 
Ito rule 



dWfrdWL, = [SijSi.j, - S ir S Vj ]dt 
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Following a standard procedure JTJ] it is possible to obtain the Fokker-Planck equation (FPE) that is mathematically 
identical to the SDE (0). The result is 
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where pt(e) is the density of probability of finding a particular distribution e\, . . . , e^r of internal energies in each 
particle. 

It is noteworthing that if one assumes the detailed balance condition a 2 = 2k, and u>(r) = oj 2 (r), then the FPE 
becomes 
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d tPt (e) = L^p t (e) 

where the heat conduction operator is given by 

y 

The FPE (|) admits the following unique equilibrium solution 

Peq(e) = iexp|^s(e 4 )|p(^£fc) 

where Z is the normalizing factor and VQ2 k efc) is a function of the total internal energy and will be discussed in the 
next section. This is the Einstein's formula for energy fluctuations in the presence of energy conservation. Note that 
at equilibrium the intensive parameters T~ Y are equal on average [fl5[ 

57) =/ deA[eKpS(e)]7>(5>) 
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On the other hand, if one computes with the method of Lagrange multipliers the most probable distribution efj, . . . , ejv 
subject to the conservation of internal energy one also obtains T(?i) = T(ej) |l5| . 



III. VISCOUS HEATING 



When the particles have certain velocities v^, the DPD model prescribes the following equations of motion 
dY; — \f;dt 
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where is a conservative force that derives from a potential Vfcj), v%j — Vj — Vj is the relative velocity, and the 
functions B(r),B(r) provide the range of the forces. In the original DPD algorithm the friction coefficient 7^ and 
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the noise amplitude Oij where assumed to be identical for all pairs, this is, 7^ = 7 and er^ = a. We need to be more 
general now due to the temperature variations. The increments of the Wienner process are symmetric under particle 
interchange, dW^ = dW^, in order to ensure momentum conservation j?J. They satisfy the Ito mnemotechnical rule 



(10) 



From the equation of motion (g) it is possible to evaluate an infinitesimal increment of the total mechanical energy 
defined as E mec — | Y^[j ^fajO + Et \ rav 'i- ^ ne P rmie m the summatory excludes the terms i = j. Using Ito calculus 



with the rules (|10[) we obtain 
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We make the assumption that all the dissipated energy (|llj) is invested in rising the total internal energy E — — e, 
in such a way that the total energy E + e$ is exactly conserved. This can be achieved simply by considering the 
following SDE for the internal energy 
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Now, if we assume that both processes of conduction and viscous heating are operating then the proposed equations 
of motion are given by 



dri = Vjdt 
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We postulate that dWZ and rfW^ are uncorrelated. 

These stochastic differential equations have associated an equivalent Fokker-Planck equation that governs the evo- 
lution of the distribution function p t of the all the variables of the system. By assuming 



ui 2 (r) = ui(r) 

a 2 = 2k 
B 2 (r) = B(r) 

and after some straightforward although lengthy algebra the following FPE results 

d tPt =L FF p t =[L c + L VH +L HC ] (h 
where L c is the usual Liouville operator of the conservative system 
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The heat conduction operator L HC is given by (o) and the viscous heating operator L VH is given by 
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where the operator Ly is defined by 



Lij — ejj • 
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It is most important to realize that the Fokker-Planck operator in (15) satisfies L FP f(E) — and L FP g(P) = 
for any arbitrary functions /, g of the total energy and total momentum. This crucial property is a reflection of the 
conservation of energy and momentum of the SDE ( |l3| ) . 

As in the previous section, the requirement that the FPE has a prescribed equilibrium solution imposes certain 
constrains on the different functions and parameters of the model. We postulate that the equilibrium distribution 
function is given by 



Peq 



(r, v,e) = - exp{J2 s{e t )}P(E tot (r, v, e), P) 
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where P(E, P) is a function of the total energy and momentum to be determined below. It is easily shown that if 
<7y = cr, a constant for all particles, and 
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then the equilibrium distribution function ( |l9| ) is the solution of the FPE (|5|) . This non-trivial result ( pp| ) constitutes 
along with ( |l4| ) the detailed balance conditions for the model. The result ( (20| ) is to be compared with the detailed 
balance condition of the algorithm without internal energy, 7 = ma 2 /2T. 

We discuss now the meaning of P(E,P) by considering the time-dependent probability density P(E,P) that the 
system has a given value of the dynamical invariants E,P. It is given by 



V(E,P-t) = J dzS(E - E(z))S(P -^mv fc )ft( 2 
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where z is a shorthand for all the variables of the system and E{z) — \ Ylij V( r ij) + J2i \ mv 1 + Si e i- The time 
derivative of (|Tj) is zero, d{P = 0, as can be seen after using the Fokker-Planck equation for dtpt(z) and an integration 
by parts. This means that the equations of motion do not modify the distribution of dynamical invariants, which is, 
therefore, determined by the initial distribution function po(z). By using the equilibrium distribution function p eq {z) 
in (2l]) one obtains P(E,P) = V(E,P) and the function P(E,P) in ( |l9| ) is determined by the initial distribution 
of dynamical invariants. A similar argument can be given in classical mechanics where instead of a Fokker-Planck 
operator a Liouville operator applies Jig] , In general, we will assume that P(E,P) = P(E)S(P) because the center 
of mass of the total system will be selected at rest at the initial time. 

A final comment on the meaning of the temperature is in order. The rate of change of mechanical energy is inferred 
from equation (Jll|), 
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The first contribution is the energy dissipated by friction whereas the second contribution is the energy provided by 
the noise. At equilibrium both terms are equal and there is no net variation of the mechanical energy, on average. 
These terms can be computed explicitly with the distribution function ( |ig| ) by using the invariance of ( |l9| ) under 
rotations (see appendix in Ref. p7|). It turns out then that the following relationship is satisfied (for any particle k) 



(23) 



where D is the physical dimension of space. This is a precise statement about the relationship between the temperature 
of the particle and its kinetic energy, at equilibrium. It constitutes the equivalent to the equipartition theorem to 
which it reduces if all the particles have a constant temperature and consistently the canonical distribution is assumed. 

In summary, we have presented an extension of DPD that includes an internal energy variable. The mechanical 
energy dissipated by friction is transformed into internal energy. In addition, heat conduction occurs. The macroscopic 
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behavior of the system will be hydrodynamic in the momentum and energy. It is possible to construct a kinetic theory 
along the lines of Ref. JhJ for the DPD model of this letter and obtain balance equations for the momentum and 
energy. The energy now is a slow independent variable (because it is conserved) and it is possible to have thermal 
gradients. This opens up the possibility of study a large variety of problems in which thermal processes in complex 
fluids are of concern. 

Discussion with Marisol Ripoll during the elaboration of this work are gratefully acknowledged. This work has been 
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